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Abstract 

We study the evaporation of stars from globular clusters using the simplified Chan- 
drasekhar model [S. Chandrasekhar, Astrophys. J. 97, 263 (1943)]. This is an ana- 
lytically tractable model giving reasonable agreement with more sophisticated models 
that require complicated numerical integrations. In the Chandrasekhar model: (i) 
the stellar system is assumed to be infinite and homogeneous (ii) the evolution of 
the velocity distribution of stars f{v,t) is governed by a Fokker-Planck equation, the 
so-called Kramers-Cliandrasekhar equation (iii) the velocities \v\ that are above a 
threshold value R > (escape velocity) are not counted in the statistical distribution 
of the system. In fact, high velocity stars leave the system, due to free evaporation or 
to the attraction of a neighboring galaxy (tidal effects). Accordingly, the total mass 
and energy of the system decrease in time. If the star dynamics is described by the 
Kramers-Chandrasckhar equation, the mass decreases to zero exponentially rapidly. 
Our goal is to obtain non-perturbative analytical results that complement the semi- 
nal studies of Chandrasekhar, Michie and King valid for large times t +oo and 
large escape velocities R — > +oo. In particular, we obtain an exact semi-explicit solu- 
tion of the Kramers-Chandrasekhar equation with the absorbing boundary condition 
f{R,t) = 0. We use it to obtain an explicit expression of the mass loss at any time 
t when R — > +oo. We also derive an exact integral equation giving the exponential 
evaporation rate X{R), and the corresponding eigenfunetion fx{v), when t — > +oo for 
any sufficiently large value of the escape velocity R. For R — > +oo, we obtain an ex- 
plicit expression of the evaporation rate that refines the Chandrasekhar results. More 
generally, our results can have applications in other contexts where the Kramers equa- 
tion applies, like the classical diffusion of particles over a barrier of potential (Kramers 
problem). 



1 Introduction 



In a seminal paper, Chandrasekhar [T] developed a Brownian theory of stellar dynamics in 
order to determine the rate of escape of stars from globular clusters. Small groups of stars 
tend to approach a statistical equilibrium state (described by the Boltzmann distribution) 
as a result of stellar encounters. However, high energy stars are not bound to the system 
and escape to infinity. For an isolated system, the average escape velocity for all stars 
in the cluster is fixed by the virial theorem according to the equation R = 2vrms where 
Vrms = {v^)^/"^ is the root-mean-square velocity (RMS) [2]. If the system, e.g. a globular 
cluster, is submitted to tidal forces from a neighboring galaxy, the escape velocity can 
be smaller. Therefore, stars clusters tend to slowly evaporate. This evaporation was first 
studied by Ambartsumian [3] and Spitzer [4] using phenomenological arguments. They 
estimated the evaporation time by removing a fraction 7 = 7.38 10"'^ of stars every re- 
laxation time, where 7 is the fraction of particles in a Maxwellian distribution that have 
speeds exceeding twice the RMS velocity. In a more precise treatment, Chandrasekhar 
[1] described the "collisional" evolution of a stellar system by a Fokker-Planck equation 
(the nowadays called Kramers equation) involving a diffusion term in velocity space mod- 
eling the erratic motion of the stars and a friction term that appears to be necessary to 
drive the system towards the Boltzmann distribution predicted by statistical mechanics 
(fluctuation-dissipation theorem). The diffusion coefficient and the "dynamical friction", 
satisfying the Einstein relation, were independently justified from kinetic theory by ex- 
plicitly calculating the first and second moments of the velocity increment suffered by a 
star during a succession of binary encounters. In order to account for evaporation, Chan- 
drasekhar imposed as a boundary condition that the distribution function f{v,t) vanishes 
when the star velocity reaches a maximum value \v\ = R. He then reduced the problem 
to the study of an eigenvalue equation in a bounded domain of velocities \v\ < R. The 
fundamental eigenvalue gives the exponential evaporation rate of the stars from the clus- 
ter for t —7- +00 and the associated eigenfunction gives the quasi-stationary distribution 
function of the system. This distribution function is close to the Boltzmann distribution 
for low velocities and tends to zero at the escape velocity. Chandrasekhar solved the 
eigenvalue problem by transforming the Kramers equation into a Schrodinger equation 
(with imaginary time) for a quantum oscillator in a box and by expanding the solutions of 
that equation in the form of a series. He obtained (semi-explicit) analytical results in the 
R —7- +00 limit or, equivalently, for a small evaporation rate X{R) — )• 0. In his first treat- 
ment [1], he assumed that the diffusion coefficient is constant and in a more exact theory 
[5], he took into account the dependence of the diffusion coefficient with the velocity. His 
work was followed by Spitzer & Harm who determined the escape rate (eigenvalue) 
and the quasi-stationary distribution function (eigenfunction) numerically for any value 
of the escape velocity R. Then, Michie [7] and King [8] obtained for R — ?• +00 a simple 
analytical expression of the quasi-stationary distribution function in the form of a lowered 
isothermal distribution which vanishes at the escape velocity. This leads to the so-called 
Michie-King model that is asymptotically valid in the limit R — >• +00. 
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The Chandrasekhar model described previously is based on simplifying assumptions. 
It is first assumed that the system is spatially homogeneous and infinite while globular 
clusters are highly inhomogeneous and limited in space. On the other hand, the coUisional 
evolution of the system is modeled by the Kramers equation while a more relevant equa- 
tion is the gravitational Landau equation that is the standard kinetic equation of stellar 
dynamics [H [9l [TU[ The Kramers equation corresponds to a canonical description 

in which the system is assumed to be in contact with a thermal bath from which it can 
extract energy so that the temperature T is fixed. Alternatively, the Landau equation 
corresponds to a microcanonical description in which the system is assumed to be isolated 
so that the energy E is conserved. Since globular clusters are isolated Hamiltonian sys- 
tems (up to the slow evaporation process), the microcanonical description appears to be 
more relevant. Therefore, when we take into account spatial inhomogeneity and model 
the encounters in a self-consistent way, the proper model to consider is formed by the 
gravitational Landau equation coupled to the Poisson equation. In order to go beyond 
the limitations of the Chandrasekhar model and obtain more accurate rates of escape, the 
astrophysicists have performed numerical simulations of stellar systems. Different types 
of simulations have been performed. They solved (i) the A^-body Hamiltonian problem 
associated to the Newton equations [12] (ii) the hydrodynamic moments of the Landau 
equation [131 [T^ (iii) the A-body problem where the effect of encounters is modeled by 
Monte Carlo methods [151 [TBI I17j . and (iv) the orbit-averaged- Fokker-Planck equation 
|18j . These methods are reviewed in the books of Spitzer [9] and Binney &: Tremaine [2] 
and in the reviews [19^ 120] . In these numerical works, the spatial inhomogeneity of the 
cluster is properly taken into account. 

These simulations have led to the following scenario for the evolution of globular clus- 
ters [21 |9]. In a first regime, a self-gravitating system initially out-of- mechanical equi- 
librium undergoes a process of violent coUisionless relaxation towards a virialized stat^. 
In this regime, the dynamical evolution of the cluster is described by the Vlasov-Poisson 
system and the phenomenology of violent relaxation has been described by Henon |21j . 
King [22] and Lynden-Bell [23]. Numerical simulations that start from cold and clumpy 
initial conditions generate a quasi stationary state (QSS) that fits the de Vaucouleurs 
R^/^ law quite well [24]. The inner core is almost isothermal (as predicted by Lynden-Bell 
[23]) while the velocity distribution in the envelope is radially anisotropic and the density 
profile decreases like [251126]. One success of Lynden-Bell's statistical theory of violent 
relaxation is to explain the isothermal core without recourse to "collisions". By contrast, 
the structure of the halo cannot be explained by Lynden-Bell's theory as it is a result of 
an incomplete relaxation. On longer timescales, encounters between stars must be taken 
into account and the dynamical evolution of the cluster is governed by the Vlasov-Landau- 
Poisson system which is the standard model of stellar dynamics. This collisional regime 
is appropriate to understand the actual structure of globular clusters. In this regime, the 

^This form of relaxation is appropriate to account for the actual structure of elliptical galaxies whose 
dynamics is encounterless for the timescales of interest [2]. 
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system passes through a succession of quasi stationary states (QSS) that are steady states 
of the Vlasov equation slowly evolving in time due to the cumulative effect of encounters. 
The first stage of the collisional evolution is driven by evaporation. Due to a series of 
weak encounters, the energy of a star can gradually increase until it reaches the local 
escape energy; in that case, the star leaves the systerrH. Numerical simulations [9] show 
that during this regime the system reaches a quasi-stationary state that slowly evolves in 
amplitude due to evaporation as the system loses mass and energy. This quasi stationary 
distribution function (DF) is close to the Michie-King model. The system has a core-halo 
structure. The core is isothermal while the stars in the outer halo move in predominantly 
radial orbits. Therefore, the distribution function in the halo is anisotropic. The density 
follows the isothermal law p ~ in the central region (with a core of almost uniform 
density) and decreases like p ^ r^''/^ in the halo [2]. Due to evaporation, the halo expands 
while the core shrinks as required by energy conservation. At some point of the evolution, 
when the energy passes below a critical value (or when the density contrast becomes suf- 
ficiently high), the system undergoes an instability related to the Antonov |28| instability 
and the gravothermal catastrophe takes place [29]. This instability is due to the nega- 
tive specific heat of the inner system that evolves by losing energy and thereby growing 
hotter (see reviews in [10l[30]). This leads to core collapse [2]. Mathematically speaking, 
core collapse would generate a finite time singularity: if the evolution is modeled by the 
orbit-averaged-Fokker-Planck equation, Cohn |18] finds that the collapse is self-similar, 
that the central density becomes infinite in a finite time and that the density behaves like 
p ~ 7--2.23 ^YiQ evolution is modeled by the Landau-Poisson system, it is argued in [31] 
that the density behaves like p oc in the final stage of the collapse). In reality, if we 
come back to the A^-body system, core collapse is arrested by the formation of binary 
stars. These binaries can release sufficient energy to stop the collapse j32j and even drive 
a re-expansion of the cluster in a post-collapse regime [33]. Then, in principle, a series of 
gravothermal oscillations should follow [3l]. In practice, the processes of evaporation and 
core collapse take place simultaneously so that it is difficult to isolate the effect of any 
single process in the evolution of a globular cluster. 

Concerning the evaporation process, the Princeton code was the first code to yield 
reliable evaporation rates [35] giving tevap = —N{dN/dt)~^ ~ 300trh for isolated clus- 
ter^. These results are not very different from those obtained with the spatially ho- 
mogeneous Kramers-Chandrasekhar equation. Our goal in this paper is not to make a 
realistic modeling of stellar systems but rather to consider simple models of evaporation 
that are analytically tractable. Therefore, we shall use the Chandrasekhar model which 
yields a reasonable description of the evaporation process in globular clusters and which 
can be studied analytically. Chandrasekhar solved the problem perturbatively: he first 
considered the long time limit t — )■ +oo so that the distribution function fniv, t) is domi- 

^There can also be a process of ejection [27j in which a single close encounter produces a velocity change 
that is sufficient to eject the star out of the cluster. However, it can be shown that this process is less 
efficient than evaporation. 

^Tidal forces from the Galaxy can increase the evaporation rate [36| . 
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nated by the contribution of the fundamental eigenmode fR{v)e~^^^^^^^ and then took the 
hmit R — >• +00 to obtain an approximate expression of the quasi-stationary distribution 
(fundamental eigenfunction) and escape rate X{R) (fundamental eigenvalue). In 
this paper, we shall reconsider the Chandrasekhar problem on a new angle which allows 
to obtain non-perturbative results. In particular, we find an exact semi-explicit solution of 
the Kramers equation with boundary condition f{v,t) = when \v\ = R. This solution 
f{v,t) depends on the remaining mass in the cluster Mq^I) which satisfies an autonomous 
equation. We use this general formula to obtain: (i) the mass MQ{t) for any fixed time t 
in the limit R — ?• +00, (ii) an exact integral equation for the eigenvalue \{R) of the fun- 
damental mode (evaporation rate) valid for any sufficiently large i?, (iii) an exact explicit 
expression of the fundamental eigenfunction valid for any sufficiently large ii, and (iv) an 
explicit asymptotic expression of the evaporation rate when R — )■ -|-oo. Therefore, our 
approach complements Chandrasekhar's original work and offers new perspectives. Our 
main results are, however, restricted to the Kramers equation, i.e. a Fokker-Planck equa- 
tion with constant diffusion coefficient and quadratic potential (linear friction). A different 
approach that goes beyond these limitations (but which is restricted to the asymptotic 
limits t —7- +00 and R — )■ +00) is developed in Appendix [Dl 

Let us finally note that our approach is not limited to the astrophysical problem 
mentioned above but that it can have applications in different area. First of all, Chan- 
drasekhar's study of the rate of escape of stars from globular clusters is closely connected 
to the classical Kramers [S^ problem for the escape rate of a Brownian particle across 
a potential barrier that has many applications in physics and chemistry (surprisingly, 
Chandrasekhar [1] did not mention this connection). In that case, the problem is usually 
formulated in d = 1 dimension. On the other hand, Chandrasekhar's procedure has been 
used in the context of planet formation [38] in order to determine the rate of escape of dust 
from large-scale vortices (assumed to be present in the solar nebula) due to turbulence. 
In that case, the problem is two-dimensional. In view of the fundamental nature of the 
mathematical problem, it seems relevant to develop our formalism in arbitrary dimension 
of space d in order to cover a wide range of possible applications. 
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2 Setting of the problem and statement of the results 



2.1 Kinetic models on a bounded velocity domain 

Basically, the evolution of the distribution function f{x, v, t) of a stellar system is described 
by the Vlasov-Landau-Poisson equation [21 [9l [ini IH] : 

^ + w • Vxf - V^> • Vvf = Vv- j K{v - -y*) (/(x, -y*)V^,/(x, v) - f{x, v)Vi,f{x, v*)) dv^ 
A$ = iirG j f dv, 

f{x, v,t = 0) = fo{x, v) > 0, (x, v) eM.^ X M^, 

(2.1) 

where is the gravitational potential and K{u) is the following 3x3 matrix 

K{u),, = ^H^fc^, (2.2) 

where A = 27rG^mln A is a constant {G is the gravity constant, m the mass of a star and 
In A the Coulomb logarithm). These equations must be complemented by the boundary 
condition f{x,v,t) = if e = + ^{x,t) > which expresses the fact that stars with 
positive energy are lost by the system. Indeed, they are unbound and free to escape to 
infinity. This is the reason for the evaporation of the star cluster. The Vlasov-Landau- 
Poisson system (12.10 is the standard model of stellar dynamics. From the Landau equation, 
the relaxation time due to two-body encounters can be estimated by tji ~ {N/\nN)t£) 
where is the dynamical time and the number of stars in the system [2]. For large 
groups of stars like elliptical galaxies (A^ ~ 10^^, to ~ 10^ years, age ~ 10^ years), the 
relaxation time is much larger than the age of the universe by several orders of magni- 
tude so that star encounters are completely negligible. In that case, their evolution is 
described by the Vlasov-Poisson system, i.e. by (12. ip with the r.h.s. taken equal to zero. 
For smaller groups of stars like globular clusters (A^ ~ 10^, to ~ 10^ years, age ~ 10^*^ 
years) whose ages are of the same order as the relaxation time, the encounters must be 
taken into account. The study of the full Vlasov-Landau-Poisson equation is extremely 
complicated because it involves several processes: (i) violent relaxation in the collision- 
less regime (ii) gravothermal catastrophe in the collisional regime, and (iii) evaporation. 
Furthermore, the coupling between the Landau equation and the Poisson equation, and 
the fact that the distribution function depends on seven variables (x, v, t) in the general 
case, make these equations untractable without further assumptions. Some simplification 
can be obtained by averaging over the orbits of the stars thereby obtaining the orbit- 
averaged- Fokker-Planck equation [21 [9]. For spherical systems, this leads to an equation 
for the distribution function / = /(e, t) that depends only on the energy e = f ^/2 + and 
time t. Still, the theoretical study of the orbit-averaged-Fokker-Planck equation remains 
very complicated. In order to distinguish the contribution of each process occurring in 
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the evolution of a stehar system and to study specificahy the evaporation process in a 
very simple setting (which is our motivation here) we shall make additional simplifying 
assumptions. First, we shall disregard the spatial structure of the system and assume 
that the medium under consideration is infinite and homogeneous. If we implement this 
approximation naively, solving the spatially homogeneous Landau equation for f{v, t), it is 
found 09] that the r.m.s. velocity decreases due to evaporation while in reality (i.e. when 
spatial inhomogeneity is retained and the Landau equation is coupled to the Poisson equa- 
tion) the contraction of the core causes the r.m.s. velocity to increase (as potential energy 
is converted into kinetic energy faster than the escaping stars carry energy away). This 
approximation leads therefore to unphysical results. This problem was solved by King [39] 
by adding artificially in the spatially homogeneous Landau equation an additional outward 
flux in velocity space. In that case, the decrease in kinetic energy due to evaporation is 
compensated by the increase in kinetic energy due to contraction. Another solution, that 
was proposed earlier by Chandrasekhar [1] , is to assume that the star under consideration 
has encounters with a separate group of stars having a fixed (usually assumed Maxwellian) 
velocity distribution. In that case, the star is able to extract energy from a reservoir im- 
posing its temperature (thermal bath). Therefore, Chandrasekhar models the evolution 



where -D(|?;|) is some given nonnegative diffusion matrix. The Fokker-Planck equation 
(j2.3p can be derived from the Landau equation (|2.ip by making the so-called "thermal 
bath approximation", i.e. by replacing the function /(t;*) in (12. Ih by the Maxwell dis- 
tribution f{v) = /o(27rT)~^/^exp(— 'y^/2T) with inverse temperature /3 = 1/T (here, the 
mass of stars has been included in the temperature). This procedure transforms an in- 
tegrodifferential (Landau) kinetic equation into a differential (Fokker-Planck) equation 

* Chandrasekhar [I] did not give a precise justification for using a differential Fokker-Planck equation 
instead of an integro-differential equation like the Landau equation. The Fokker-Planck equation considered 
by Chandrasekhar describes the evolution of a "test star" in a bath of "field stars" at statistical equilibrium 
with fixed temperature (canonical description). By contrast, the Landau equation describes the evolution 
of the system as a whole and conserves the energy (microcanonical description). Apparently, the Landau 
equation was not well-known by the astrophysical community at that time (see discussion in [11]). Here, 
we shall assume that the existing contraction of the cluster has the effect to "heat up" the stars and 
balance their "cooling" due to evaporation. As a result of these two antagonistic effects, we can assume 
that the temperature (velocity dispersion) of the stars remains fixed. Therefore, everything happens as if 
the system were in contact with a thermal bath, justifying the Chandrasekhar assumptions. Accordingly, 
when we disregard the spatial structure of the system (but take it into account indirectly as explained 
above), it is more appropriate to model the dynamics of the star cluster by the Fokker-Planck equation 
(canonical) rather than by the Landau equation (microcanonical). However, if we take into account the 
spatial structure of the system, the best kinetic description is the Landau equation coupled to the Poisson 
equation. 




-^iv,t) = Q''/{f){v) = V ■ [D{\v\) (Vfiv) + I3f{v)v)] ,veBn, 

f{v,t = 0) = fo{v)>0, v^Br 
f{v,t) = 0, ii\v\=R, 



(2.3) 
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and allows an explicit computation of the diffusion coefficient -D(|u|) (see [TT] and refer- 
ences therein). The Fokker-Planck equation (|2.3p can also be derived from a stochastic 
Langevin equation with a linear friction (Ornstein-Uhlenbeck process) and is called the 
Kramers equation [37]. Since it was derived independently by Chandrasekhar [1] in the 
astrophysical context, we will call it here the Kramers-Chandrasekhar equation. Finally, 
the last condition in (|2.3p expresses the fact that a star with velocity above a certain 
threshold R escapes and is therefore lost by the system. This threshold can be estimated 
by the following classical argument [2]. The escape speed at x is given by = —2^{x), 
corresponding to e = 0. The mean square escape speed in a system whose density is p{x) 
is therefore (fg) = / p{x)v1dx/ J p{x)dx = — (2/M) J p{x)^{x) dx = —4W/M where M 
is the total mass and W the potential energy. According to the virial theorem —W = 2K 
where K = {1/2)M {v'^) is the kinetic energy. Hence {vD = 4(v^). For a Maxwellian 
distribution with temperature T, we obtain R = 2\/3T. However this is just an estimate. 
Furthermore, the escape velocity can be smaller if the system (globular cluster) is sub- 
ject to the tide of a nearby galaxy. Therefore, for sake of generality, we shall consider R 
arbitrary. The Fokker-Planck equation ()2.3p with proper boundary condition f{R,t) = 
was used by Chandrasekhar [1] and others [6l [71 [8] to determine the rate of escape of 
stars from globular clusters. We shall call it the Chandrasekhar model. This is the model 
studied in the present paper. For sake of generality, we shall consider these equations in 
d dimensions. 

2.2 Semi-explicit solutions on a bounded velocity domain 

In this section, we focus on the Fokker-Planck equation (j2.3p with D = 1 and boundary 
condition f{R,t) = 0. We prove that a semi-explicit solution of this kinetic equation can 
be given for this model. More precisely, we derive an autonomous relation satisfied by 
Mo(t), the mass remaining in the cluster at time t, and then give an explicit expression 
of the solution f{v,t) in terms of Mo(f) and of the initial data /o(f) only. This result 
will be used to determine the exact (i.e. non perturbative) rate of escape of stars in the 
Chandrasekhar model (see section 2.3). 

Note that the difficulty here comes from the boundedness character of the velocity 
domain and the presence of the term vf in the drift-diffusion model ()2.3p . which prevent a 
direct use of Fourier analysis. Our strategy is to first transform Eq. ()2.3p into equivalent 
equations on the whole domain where the boundary condition naturally disappears and 
is replaced by a source term (see Appendix E]). The resulting equation makes sense in 
the space of distributions and one can use Fourier techniques in this space to work out a 
semi-explicit solution of the problem. 

Proposition 2.1 (Semi-explicit solution of (j2.3p ) Let /o be a smooth and isotropic 
initial data supported inside Bn. The solution to ()2.3p when D = 1, with initial data /o, 
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is given by 

+7T^^ / ds -r-— M'^{s) / da exp ' ' ^ ' ^ 



(2.4) 

where §f zs i/ie unzi sphere with measure do" is i/ie surface element of this unit 
sphere, and: 

A{s, t) = ^ (1 - exp(-2/3(i - s))) , i?(s, t) = exp(-/3(i - s)). 
T/ie total mass Mo{t) of f, defined by 

Mo{t)= [ f{v,t)dv, (2.5) 

J Br 

is determined from the boundary condition f{Ruj,t) = 0, VtJ € Sf leading to the au- 
tonomous equation 

A(0, t) j^^ /o(..)exp ^ j d.. 

+ f dsAis,t)-^I^M',is) [ da exp f-^^^^/l^) = 0. 
Jo Jaesf V 4A{s,t) J 

(2.6) 

A'^oie i/iai this last relation does not depend on u and in particular can be averaged over 
uj G Finally Mq(s) is the derivative of Mq{s) at time s. 

The proof of this result is given in Appendix |X1 Note that the solution (j2.4|) is semi- 
explicit in the sense that it involves the quantity Mq{s) which depends on the solution itself. 
However, the only knowledge of the time evolution of this mass allows the determination 
of the whole solution /. The exact equation satisfied by Mq{s) is given by (|2.6p and, as 
it stands, seems to be too complicated for practical use. Nevertheless, this equation can 
be simplified in the asymptotic limit R — ?• +00 for fixed time t. This is the subject of the 
following proposition 

Proposition 2.2 (Approximate mass law for large R) Let f be a smooth enough so- 
lution of ()2.3p with D = 1 and isotropic initial data /q supported on Br. Let Mo(t) be the 
total mass at time t given by ()2.5p . Then, for any given time t > 

Y^7r(l -exp(-2/3t)) V ^ / ^2 7) 

^ — ^ /3(exp(-/3i)r -i?)^^ 



d-l 



r 2 fo(r)exp ; ^ dr, 

^ ' \ 2(1 - exp(-2/3i)) i 
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as R goes to +00. Furthermore, 



If fo does not depend on R and is supported on [0, Rq] with a non-vanishing left 
derivative /Q(i?o) at Rq, then 

M'^{t) 2— |Sf| ^ R^exp ( '^Pt] smh2(/3t) 

^ V7r/5(1 - exp(-2/3t)) ° V 2 7 

xexp 



/? (exp(-/3t)iio - -R)^ 
2(1 - exp(-2/3t)) 




as R ^ +00. 

// /o depends on R with a non vanishing left derivative f'^iR) at R, then 



2d+l , 



Kit) - ^^7p=2L=«'-'exp { 'LApt 



Vvr/3(l-exp(-2/3t))" 



.exp f-^(^--P^^f\)%- 1 r,iR). 



(2.9) 



(1 - exp(-/3i))2 ^ ^ \ 2(1 - exp(-2/3t)) 



as R ^ +00. 



The proof of this result is given in Appendix |Bl 

Remark: In practice, for sufficiently large R, this expression is 'valid' for t <^ R^ so 
that MQ{t) <C 1. This clearly shows that the order of the limits R — )• +00 and t — t- +00 is 
not interchangeable. Usually, most works [H O [H [71 [8] consider the limit t — ?• +00, then 
the limit R — )■ +00. By contrast, the above expressions are valid for R — )• +00 at any fixed 
time t. 



2.3 Exact rate of escape for ( 12. 3 p 

In this section, we focus on the Fokker-Planck equation (j2.3p with D = 1. It is well known 
that the long time behavior t — )• +00 of the solution to (|2.3p can be described from the 
knowledge of the first eigenvalue (the largest nonzero eigenvalue) of the linear operator 
Qjl^ constrained with the vanishing Dirichlet boundary condition. This is essentially a 
consequence of the self-adjointness of this operator in the space {exp{\v\^ /2)dv) . How- 
ever, the analytical determination of this eigenvalue is a difficult task. To our knowledge, 
the first work on this subject goes back to Chandrasekhar [1] but the expression obtained 
for the fundamental eigenvalue is only an approximation in the limit of large R and is 
given in terms of a (not explicitly summable) series. Here, we propose another strategy 
which is based on the semi-explicit solution (j2.4p and derive an exact autonomous relation 
satisfied by this first eigenvalue for any sufficiently large R. Then, as a consequence, we 
recover the Chandrasekhar result (in a more explicit form) by taking the leading term in 
our relation when R is large. Here is the statement 
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Proposition 2.3 (Exact rate of escape for (|2.3p ) Let X{R) be the largest non zero 
eigenvalue of the linear operator Qfp given by (12. Sp with D = 1, in the space of isotropic 
functions of (exp(|t>p/2)(iu) , vanishing at the boundary. Then 

i) X = A(i?) is negative and, for sufficiently large R so that X{R) + 2/3 > 0, it satisfies the 
following nonlinear relation 

G{Ruj,0) + ^ [ {G{Ruj,u) -G{Ruj,0))u^~^du = 0, V cj e Sf, (2.10) 
P Jo 

where 

for all V ^ and u G [0, 1[. Note that relation (|2.10p is independent of oo and can 

therefore be averaged over w G §f . 

a) The corresponding eigenf unction is exactly given by 

fx{v) = Mx {civ, 0) + ^ [G{v, u) - G{v, 0)] du] , (2.12) 



PJo 

with Mx = !br /A(^)f^'y- 

Hi) The eigenvalue \{R) has the following asymptotic behavior 



2/3 (pR^Y'^ (-m\ 



r(i) 

as R goes to +oo. 

The proof of this result is given in Appendix [Cl 

Remark 1: The exphcit asymptotic behavior p.lSp was not given by Chandrasekhar 
[1] who obtained the asymptotic expression of the eigenvalue in the form of a series. In 
Appendix |Dl we obtain the asymptotic behavior (I2.13P by a different method which can 
be extended to the case of a Fokker-Planck equation with an arbitrary diffusion coefficient 
-D(|t;|) and potential C/(|t;|). However, the method in Appendix [D] is formal and, as it 
stands, cannot be considered as complete mathematical proof of the result, unlike the 
proof given in Appendix O 

Remark 2: The expression of the function G{v, u) can be simplified as shown in Ap- 
pendix |e1 

3 Numerical results and discussion 

We have performed numerical simulations in order to illustrate our theoretical results. We 
have taken d = 3 (appropriate to stellar systems) and D = P = 1. 
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Figure 2: Normalized velocity profile f{v)/f{0) at large times (eigenfunction) for different 
values of R. The slope at the escape velocity v = R decreases as R increases leading 
to slower mass loss. For R — t- +oo, the distribution function tends to the Maxwellian 
(bullets). 
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Figure 3: Decay of the mass Mo(t) contained in the cluster as a function of time (logarith- 
mic scale). For large times, the decay is exponential Mo(t) ~ e~^^^^'^^^ leading to straight 
lines with slope X{R). The exponential decay rate |A(i?)| decreases as R increases. 




' 



0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 

R 



Figure 4: Mass decay rate |A(i?)| (eigenvalue) for different values of R. The decay rate 
tends to +00 for — )• and is equivalent to the asymptotic expression ()2.13p for R — )• +00. 
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In Fig. [H we show the velocity distribution f{v, t) at different times. We have adopted 
the value R = 2\/3 of the escape velocity corresponding to the estimate deduced from the 
virial theorem (see Introduction). The initial distribution is fo{v) = e~^'" for v < 1 
and /o(f) = for ?; > 1. Due to evaporation, the distribution function decreases and 
tends to zero for t — t- +oo. In fact, its large time behavior is of the form f{v,t) ~ 
g-|A(_R,)|ty^^^^_ If -vvg rescale the distribution f{v,t) by its central value f{0,t), then the 
normalized velocity profile f{v,t)/f{0,t) tends to the eigenfunction fx{v)/fx{0). This 
eigenfunction is represented in Fig. [2] for different values of R. We see that the slope 
of the distribution at the escape velocity R decreases as R increases implying a slower 
mass loss. In fact, for R — t- +oo, the mass loss tends to zero and the eigenfunction tends 
to the Maxwellian which is the steady state of the Kramers equation without velocity 
confinement. In Fig. [3l we plot the mass Afo(t) contained in the cluster as a function 
of time in logarithmic scale. For large times, the decay is exponential Mo(t) ~ e"'"^^^^'* 
leading to straight lines with slope A(i?). The exponential decay rate |A(i2)| decreases as 
R increases. This decay rate is plotted in Fig. |3]for different values of R. The decay rate 
tends to +oo for i? — )■ and behaves like ()2.13p for R — )■ +oo. 

4 Conclusion 

In this paper, we have obtained new analytical results for the escape of stars from glob- 
ular clusters in the framework of the Chandrasekhar model. These results can also have 
applications to other physical systems described by the Kramers equation with parabolic 
potential and absorbing boundary conditions (i.e. the classical Kramers problem). 

We have first obtained an autonomous equation for the mass loss [see (|2.6p ] and a 
semi-explicit expression of the distribution function f{v,t) [see (12.40 ] that are valid for an 
arbitrary escape velocity R and time t. We have simplified the expression of the mass loss 
in the limit R — )■ +oo for any fixed time t [see (12. 7p ]. We have also used the semi-explicit 
expression of the distribution function f{v,t) to obtain an exact integral equation for the 
fundamental eigenvalue A(i?) [see (|2.10p ] and for the fundamental eigenfunction fx{v) [see 
(j2.12p ] that are valid for any sufficiently large R. This is an interesting complement to the 
perturbative results derived by Chandrasekhar [l] that are valid in the asymptotic limit 
R — )• +00. We have obtained the explicit behavior of the fundamental eigenvalue in the 
limit R — )• +00 [see (j2.13p ] which improves upon the result of Chandrasekhar expressed 
in the form of a series. Additional asymptotic results for the fundamental eigenvalue 
and fundamental eigenfunction are given in Appendix [Dj Finally, we have illustrated our 
results with numerical simulations [see Sec. [3]. 

Of course, our approach is based on several approximations. As previously discussed, 
it assumes that the medium is infinite and spatially homogeneous and that the encoun- 
ters between stars can be described by the Kramers equation (Chandrasekhar 's model). 
Furthermore, our analytical results (except those of Appendix |D]) are valid only when the 
diffusion coefficient in the Fokker-Planck equation is constant while a more exact descrip- 



14 



tion of encounters between stars would involve a velocity dependent diffusion coefficient 
O El . It remains therefore a challenging issue to extend the mathematical theory of the 
escape of stars from globular clusters in the case of more realistic models. 

Acknowledgement. M. Lemou was supported by the French Agence Nationale de la 
Recherche, ANR JC MNEC. 

Appendix 



A Proof of Proposition 12.11 



We consider the Fokker-Planck equation (j2.3p with D = 1. We recall that the initial data 
is assumed to be isotropic (i.e. it depends on the modulus of the velocity only), which 
implies that the solution f{v,t) is isotropic at any time. 

First we claim that if / is a solution to (j2.3p on the ball Br with a vanishing boundary 
condition, then / is a solution on all in a distributional sense to 

^{v,t) = V- {Vf{v) + Pf{v)v) + M^(t)^, (A14) 



where Mq is related to / by (j2.5p and 6gd is the Dirac operator on the sphere: 

< 5§d >= / <j){v)da{v), 



for all test function (p. Note that on isotropic test functions (p, this Dirac operator reduces 
to: < 5gd,(/> >= \S%\(I){R). 

To prove (|A.14p . consider an isotropic solution / to (]2.3p which we extend to outside 
Bji, and integrate (j2.3p against an isotropic test function (p on M'^. We obtain after 
integration by parts: 



d 
It 



f{v, t)<P{\v\)dv) = [ {AcP-/3v V<P) f{v)dv +1 [ Vf-n da{v) ] <P{R). 



Taking i;^> = 1 on B^i, we also have 

M^(t) = / Vf-n da{v). 



Finally, these last two identities imply that / is a solution to ()A.14p in the sense of 
distributions. Note that the choice of isotropic test functions (p is sufficient thanks to the 
radial symmetry of /. 

Remark: We have shown that the solution of (j2.3p is solution of ()A.14p . This is 
sufficient for our purposes. However, it is easy to prove that, reciprocally, (|A.14p admits 
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an infinity of smooth solutions and that, among all of them, the solution that satisfies 
f{Roj,t) = is the solution of ([23]). 

We shall now use Fourier transform techniques to find out a semi-explicit solution to 
()A.14p . Let f{£,,t) be the Fourier transform of / in the v variable: 



fi^,t)= fiv,t)exp{-iv-^)dv, 
f{v,t)= [ f{tt)ewiiv-0 



Taking the Fourier transform in ()A.14p . we get 



dt 

where 



it t) = -ei - K ■ V/ + M'^{t)H{i), 



HiO = j^^ exp{-iv ■ da{v). 



We now use the method of characteristics to solve this equation. Let 

Easy computations yield 
OF 

-g^ic,t) = -eeMwt)F{i,t)+M'^{t)H{ieMm- 

This can be written as 



d_ 

dt 



exp (^|^(exp(2/3t) - l)j = exp (^|^(exp(2/3t) - 1)^ M'^{t)H{^eMPt))- 



Integrating over time, we get 



F(e,i) = exp (-|^(exp(2/3t) - 1)) F(e,0) 



+ i exp (^-|^(exp(2/3t) - exp(2/3s))) M'^{s)H{ieMPs))ds. 



As = F(^exp(-/3t),t), we obtain 

f{C,t) = exp - exp(-2/3t)) ) /(eexp(-/3t), 0) 



+ ^*exp - exp(-2/3(i - s))^ M^(s)F(eexp(-/3(t - s))) ds. 
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Using the notations of proposition 12.11 this can also be written as 

/(^, t) = exp {-eA{0, t)) fiCB{0, t),0) + r exp {-fAis, t)) M(,is)H{CB{s, t)) ds. 

Jo 

Now, we take the inverse Fourier transform using the identities 

InvF (exp(-^a W = pip ©*^>=P (-^) , 

and 

InvF (exp(-^e')exp {-iB^ • w)) (v) = Q^'^^'exp {^- ^""'/^ ) > 



to get the desired expression (j2.4p . 

Remark: if we integrate (j2.4p on the velocity, we find a trivial result: Mo(t) = Mo(t), 
which shows the consistency of this equation. 

B Proof of Proposition 12.21 



Let / be the solution to ()2.3p with initial data /o and let Mq be defined by (j2.5p . From 
(j2.4p and the vanishing boundary condition f{Ruj, t) = 0, w E 5^^, we get 

ri+r2 = 0, (B.15) 

where 

= \Si\A{0,t)-'/' I^JoMe^p (B.16) 

= f dsA{s,tr^l^M',{s) [ exp ("^^^J^^^^f^) da, (B.17) 

where A and B are defined in proposition 12.11 

To analyze the asymptotic behavior of the term T\ when R goes to +oo, the time t 
being fixed, we introduce a spherical system of coordinates (in M"'), write = rcj, a G Sf , 
cr • w = cos and obtain 



/' + 00 /"TT 

Ti =|Sf|^(0,t)-'^/2 / / /o(r 

JO JO 



, ii^-2ii5(0,t)rCOS^ + S(0,t)V2\ d„l , , . .^d-2,^ 

xexp( ^^^Q ^ J Crfr*^ ^dripMidf 'dO, 



where Cd is given by 



Cd= ,^ , 'i>2. (B.18) 



/ {smey-^d9 
Jo 
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In the following calculations, we shall assume d > 2 but we have checked by a specific 
calculation that the final result remains valid for d = 1. Making the change of variable 
5 = — cos 9, we obtain 



Ti= Crf|Sf|A(0,t)-'^/2exp ^^^^^^^ 
We then make the change of variable u = 2Mot) -^'''^^ ~^ ^) integral over 5, and get 



Ti = Cd2°'-2|Sf|A(0,t)-^/2 



B{0,t)R 



d-l 
2 



/o JO 

Now let R go to infinity to obtain 



2 / 



/o(r)exp( 4^(0^,) ^ 



2, (B.19) 



as R goes to +oo. 

Assume first that /o does not depend on R and is supported on [0, iio]- that case, 
we perform the change of variable: 



B(0,t)R,^ 



and get 



/ {B{0,t)r-Rf \ ^ 2A(0,t) f {B{0,t)Ro-Rf 
/o(r)exp — r 2 dr = — — exp ' 



/o V 4^(0,t) y B{0,t)R V 4A(0,t) 

^TO^ / 2A(M) \ , B{0,t)Ro ^4(M) 2\ 2^(0, t) \ — 

If we let R go to infinity in this expression, then we obtain 
/o(r) exp TTTTTT^ M ' ~ 







/2A(MV ^ iBiO,t)Ro-Rr \ ^""-'^^ 



18 



for large R. The case where /o does depend on R can be treated similarly to yield 
/'+°° „ . . / {B{0,t)r - R)' \ ^ 

(2Am_y ( {i-B{o,t))\ ,\ r"^ ^^-"'^ 

U(0,t)^y 4A(0,t) ^ j [l-i3(0,t)]2-^°^^^' 

for large We now deal with the asymptotic behavior of the term T2 given by ()B.17p 
when R goes to +00. First, we write T2 in the following form 

We then introduce spherical coordinates (in M*^) in the variable o" G and obtain 
Performing the change of variable u = ^2A{s't) ~^ cos 6) in the integral over 6, we get 

-li I u 2 exp(— li) an. 



R'^B{s,t) 



Letting R — )• +00, we have 



We now use the following change of variable in this last integral over s: 

l- B{s,t) ^ ^ //3 1-exp(-/3(t-g)) \ ^ 
27A(M) V2 1 + exp(-/3(t-s))y' 



to obtain 



d-l\ A. fKwT^^^t=m) / 2r^V'"^'''V 2r 



/3(l-cxp(-/it)) \l/2 / n 2 \ -(rf+l)/2 / r, 2 \ ^ 



This yields an equivalent of T2 



for large R. 

Combining (\B.l5\i with ()B.19P and ()B.22p leads to the result (fTT]) of proposition 
Finally, substituting (Ib:20]1 (resp. i^K2l\\ ) into ([221) yields dSS]) (resp. 1^ ). 
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C Proof of proposition 12.31 

We first prove i) in proposition 12. 31 We shall use the exact solution ()2.4p stated in propo- 
sition 12.11 Let A be the fundamental (largest non zero) eigenvalue of the Fokker-Planck 
operator (j2.3p . Even if the proof of A < is classical, we give the argument here for the 
sake of completeness. If / is an associated eigenfunction, the eigenvalue problem can be 
written as 

V- [Z?(|t;|)exp(-/5|^|V2)V(/(t;)exp(/3|^;|V2))] =A/. 
Integrating against /(u )exp(/3|'yp/2) on Br, we get 

D(|z;|)exp(-/3|7;| 2/2) | V {f{v)eMP\v\'' /2))\^ dv = A / /(t;)2exp(/3|z;| 2/2)^^;, 

Br JBr 

which implies that A < 0. We now prove relation (12.10p . Let fx{v) be an eigenfunction 
associated to A, then exp{Xt)fx{v) is the solution of (|2.3p with initial data fx- Therefore, 
using (|2.4p . the eigenvalue problem can be written as 



exp{Xt)fxiv) = mfx){v,t) + XMx I G(t;,exp(-/3(t - s))exp(As) ds, 

Jo 

where G is given by (j2.1ip . Mx = J^^ fx{v)dv and 



1 / TT f f \v-v,B{0,t)\^ 



Tr(/.)(M)=(^^^j y^/A(--)exp^- ' ]dv.. (C.23) 



This relation can be rewritten equivalently 

exp{Xt)fx{v) = Ti{fx){v,t) - MxG{v,0) + MxG{v,0)exp{Xt) 

+XMx [ [G{v, exp(-/3(t - s)) - G{v, 0)] exp(As)ds, Vt > 0. 



(C.24) 







We now analyze this relation for large time t. First, we expand the term Ti{fx){v,t) 
and obtain after easy computations: 



Ti{fx){v,t) = G{v,0) 



Mx + Pv-[ I vJxMdv, ) exp(-/3t) + \id- /3?;2)MAexp(-2/3t) 



( / \v*ffxMdvA exp(-2/3t) + ( I {v ■ v,f fx{v,)dv, ) exp(-2^t) 



+ 0(exp(-3/3t)). 



As fx is supposed to be radially symmetric in v, the order 1 contribution in exp(— /3t) 
vanishes, and we get after some rearrangements 



Tr{fx){v,t) = G{v,Q) 
+0(exp(-3/3t)). 



Mx + -id - (3v') ( Ma - ^ / \v,\'fx{v,)dv, ] exp(-2/3t) 
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This implies in particular that 



hm exp(2/30 [ri(/A)(^;, t) - MxG{v, 0)] 

t— >+oo 

is finite. Now we multiply (|C.24p by exp(2/3t) and perform the change of variable u 
exp(-/3(t - s)) in the rhs of (IC.24j) to obtain 

exp ((A + 2f3)t) h{v) = exp(2/3t) [Ti{h){v, t)- MxG{v, 0)] + M^exp ((A + 2/3)t) 



G{v,0) + 



cxp(-/3t) 



iG{v,u)-G{v,0))u^~^du 



(C.25) 



Before taking the limit t — )• +00, we first claim that 

A + 2/3 > 0, 



at least for large enough R. Indeed, we know that if i? = +00, is an eigenvalue of the 
Fokker-Planck operator (j2.3p . Therefore, the first eigenvalue A(i?) on must go to 
when R goes to infinity, and this proves the claim. Passing to the limit t — )■ +00 in ()C.25p . 
we conclude that we necessarily have 

fx{v) = Ma S^G{v, 0) + I [G{v, u) - G{v, 0)] J-^du 

which is relation (j2.12p . Finally, writing this relation at the boundary v = Rco, uj G Sf, 
and recalling that f\{Rijj) = 0, we get 



A 



G{Ru;,0) + ^l [G{Ruj,u) - G{Ruj,0)]u^ ^du = 0, 
P Jo 

which is relation ()2.10p . 

We now prove the asymptotic behavior ()2.13p for large R. As A(i?) goes to when R 
goes to infinity and 

G(i?..,0) = ^(2^/3)'^/2exp( -^^^ 

we obtain 

[(1 _ u')-d/^e^p 



-/3- 



exp(-/3^)~-|^ 
This is also equivalent to 



da 



PR^ \(TU + uj\ 
l-v? 



\{R) 



IdiR)' 



exp ( -P^ 



(C.26) 



where 



UR) = daf^ [(1 - .^)-'^/^exp - 1 



du 
u 



(C.27) 
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Now we write 
with 



da 



(1 _ u2)-d/2g^p / 



2'u{u + a ■ bj) 



du 

u 



and I~{R) = Id{R) - I^iR), S+ = {a e Sf,a ■ oj > 0}, S" = Sf - S+. First, we show 
that R~'^I^{R) is uniformly bounded in R. Indeed, we have 



l^^(^)l < 



\Sf\ 



(1 _ n2)-'^/2exp ( -PR 



du 
u 



and the change of variable 



leads to 



2 



1 + 



exp(— s ) — 1 



ds 



^(^ + 1) 



A simple computation of the derivative of the function 

„2 \d/2 



9{s) 



1-1 + 



pR^ 



exp(-s ), 



shows that it is an increasing function on [0, +oo[, for large enough R {R > a/ d/2j3). 
Therefore, it is a non- negative function on [0, +oo[, and consequently 



1 - exp(-s2) 



ds 



+ 1 



This clearly shows that R "^I^^R) is uniformly bounded for large R. 

We now focus on the behavior of 77 (R) when R goes to infinity and let 6 = —a ■ uj and 



u{5 — u) 



Then, we have 



da 



(1 - u2)-'^/2exp (pR^Ssiu)) - 1 



du 



(C.28) 



(C.29) 



We now show that the dominant part in I]^ (R) is given by the contribution at u = us, 
where us is such that 

max Ss(u) = Ssius). 

ne[0,l] 
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We first have 



us 



^ (i - - 5^) , = 1 (i - yr 



52 



and observe that 



Ssi'u) - Ss{us) 



Thus, from (10291) 



5u-l + 



2(1 -n2) (1-^1^52 



/3i?2 



(1 - n2)-'^/2g^p 



6{u- us) 
2n5(l -n2) 



exp 



/5i?2 



Sus 



du 

u 

(C.30) 

We now introduce a spherical system of coordinates (in W^) for the integration on the 
sphere and use identity 



F{-a ■ uj)da = Cd I F{- cos 9) {sin 9 f-'^dO = Cd I F{5){1 - 5'^)'^ d5.. 



where is a real function and Cd is given by (|B.18p . In the following, we assume d > 2 
but we have checked by a specific calculation that the results remain valid for d = 1. Using 
this identity in (IC.SOp . we get 



-1 



1 / /3i?2 



d5 (1 -(5^) — exp 







(1 - n2)-'^/2exp -(3R 



,2 (5 (n - us) 
2n5(l -n2) 



exp 



5us 



We now perform the following change of variables in ()C.3ip for both 6 and u 

?2\ 1/2 



du 
u 

(C.31) 



t 



u-us ( pR^ 5 ^ 



2 Us 



(^) (1-^^)^^ 



which is also equivalent to 



a + + 



,4\ 1/2 



(C.32) 
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with a = j^E? /2. After some calculations, we get from ()C.3ip 



I^{R) = 2Cd 



exp 



dr r 



2d- 3 



-1/2 



dt 



(a-r2)i/2 



(1 - n(t, r)2)-'^/2exp(-r" - t') - exp (-a) 



(l-n(t,r)2)3/2 



+ u{t,r){l-uit,r)us{r)) 

(C.33) 

We now want to analyze the asymptotic behavior of (|C.33p when a = goes to 

infinity. Thanks to the strong decreasing properties of exp(— — r^), one can expand 
the terms inside the integral in the limit a — )■ +oo. To do so, we first obtain from easy 
calculations 



1 — u{t, r) 



a L 



+ 



2a L 



O 



and 



1 - u{t,r)us{r) = -{t" + 2r2)V2 [(^2 + ^^2^/2 _ ^ 



Then, we plug these relations into (|C.33p to obtain 



I^{R) ~ 2Ca 



PR"' 



-d/2 



exp 



PR' 



+ 00 



dr 



+ 00 



dt 



xr 



2d-3(^2 ^ 2r2)-V2e^p(_^2 _ ^2) (^2 ^ 2^2)1/2 _ ^ 



(C.34) 

when a = jUR? /2 goes to infinity. To simplify expression ()C.34p we perform the cylindrical 
change of variable: t = pcos (j), r = psiiKp, with cf) G [0, vr], and get 



idiR) 



2^-^Cd 



PR' 



-d/2 



exp 1^^- J P exp(-p ) dp 



[(l + sin2(/>)i/2 + cos0] 



d-2 



(l + sin2 ,/))V2 



-sm< 



Performing the change of variable u = cos (j), we obtain 



2^~''C,, 



PR' 



-d/2 



exp 



PR 



■\ fd\ r+i [(2-n2)V2+^ 



Performing the change of variable u = \/2sint in the integral over u we get 



-1 [(2-^2)1/2^^] ^-2 

_i (2-n2)i/2 



-du = 2 



sm 



d-2 



-IT 



t + j] dt = 2'^-^ 



(sin ( 



24 



where we have set (j) = t + Tr/4 in the last integral. Using ()B.18|) . we obtain 

Therefore, I^{R) increases exponentially rapidly while -^j^(-R) increases less rapidly than 
and can therefore be neglected. Using (IC.26p . we finally obtain (I2.13p . This ends the 
proof of proposition [ 



D Asymptotic expressions of the fundamental eigenfunction 
and eigenvalue 

In this Appendix, we show that the asymptotic result (j2.13p can be directly obtained by a 
perturbative expansion of the solutions of the fundamental eigenvalue equation in powers 
of A in the limit A ^ 1 corresponding to V ^ +00 (in this Appendix, we introduce more 
physical notations and set v = \v\ G and V = R ^ W^). This method allows us to 
treat more general situations, e.g. Fokker-Planck equations with an arbitrary potential 
U{v) and an arbitrary diffusion coefficient D{v). For the Fokker-Planck equation (j2.3p 
with U{v) = /2 and D = 1, we recover ()2.13p . 

For isotropic distributions f{v,t), we consider the Fokker-Planck equation 



where D{v) and U{v) are arbitrary functions of the velocity. We consider the fundamental 
eigenmode 

f{v,t) = Ae'''g{v), (D.36) 

where ^ is a normalization constant. We can impose ^(O) = 1 without loss of generality. 
On the other hand, ^'(0) = and g{V) = 0. Substituting (lOMl) in (10^351) . we obtain the 
eigenvalue equation 

d \ d~lr^,^.^^ dg , oJU\ 



It can be integrated into 

For an unlimited range of velocities (V — )• +00), the Fokker-Planck equation ()D.35P admits 
a steady solution 

/(^,)= 1 (D.39) 
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implying that the fundamental eigenvalue is A(+oo) = 0. Now for V < +00, the eigenvalue 
X{V) < 0. However, for V ^ 1, X{V) — )• 0~ and we can formally expand the solution of 
the differential equation (1D.37P in the form 

g = go{v) + Xgi{v) + \^g2{v) + ... (D.40) 

To zeroth order, we have 



don ^ dU 



To first order in A, we get 



'¥ + ^91^ = ^TTiWr / 90{w)w<'-^dw. (D.42) 



dv ' dv v'^ ^D{v) 
With the boundary condition 50 (0) = 1, the first equation can be integrated in 

Substituting this expression in ()D.42p we obtain 

*i + ,,,|^ = ^r,-...»V-.... (D.44) 

dv ov v"- ^D(y) Jq 

that we shall solve with the boundary conditions 51 (0) = gi{0) = 0. We obtain 

9i{v) = x{v)e-^^^^^y^^'^\ (D.45) 
where is the function defined by 



X'{v) 



v'^-^D{v) Jo 



r e-^^('"^w'^-^dw, (D.46) 
Jo 



with x(0) = 0. To first order in A, the fundamental eigenfunction of the Fokker-Planck 
equation ()D.35P can be written 



9 



y) = e-/3(t^W-f^(o)) [1 + Xx(v)] . (D.47) 



Using the boundary condition ^(V^) = 0, we find that the fundamental eigenvalue is given 

by 

A(^) ~ {V ^ +00). (D.48) 

This is the main result of this Appendix. 
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If we specialize on a quadratic potential U (v) = f ^/2 (implying a linear friction like in 
(fO) '). we have 

9{v) = e-P^'/^ [1 + Xxiv)] , (D.49) 
and the function x(^) can be written more explicitly 

X'{v) = ^^.^^ e-^^w'^-'dw, (D.50) 
with x(0) = 0- For large v, we have 



1 fd\ (2\ 



d/2 



„2 



Therefore, for V — )• +oo, the eigenvalue \(y) is given by 



1^ f ^^Y'\n...-^^ 



^2) 



1 



W~-TT7^r^ my-'-. (D.52) 



which coincides with ([233]) for D = \. 
Let us give particular examples: 

(i) The King model basically, the collisional evolution of stellar systems is described 
by the gravitational Landau equation ()2.ip in d = 3 dimensions. Making a thermal bath 
approximation, i.e. replacing /(u*) by the Maxwellian distribution, and assuming that the 
distribution is spherically symmetric we obtain the Fokker-Planck equation ()2.3p with a 
diffusion coefficient (see, e.g., [TT]): 



Div) = 4 / w'e-''"T dv, (D.53) 

a,0 ' 



2 

,.2, ^ 



where X is a constant determining the relaxation time in the system. For the Fokker- 
Planck equation (j2.3p in d = 3 with the diffusion coefficient ()D.53p . the function x(^) is 
simply given by 

x\v) = ^eP'^, (DM) 



^On a strict mathematical point of view, the approach presented in this Appendix is completely formal 
since we have used a perturbative approach based on a formal asymptotic expansion. By contrast, the 
approach developed in Appendix [C] leads to an exact formula of X{V) in the case D = 1 obtained by a 
rigorous mathematical proof, and the asymptotic behavior (|2.13p is a consequence of this exact formula. 
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with x(0) = (we can appreciate the fortuitous cancelation of the integral first realized 
by King [8]). This equation is explicitly integrated and we obtain 

X{v) = ^{e^'^-l). (D.55) 

The asymptotic evaporation rate (ID.48P is then given by 

KV) = ^ i^/3e-^^. (D.56) 

el^— - 1 

The fundamental eigenfunction (|D.49P is the King solution 

p A-" 2 — p '^2 

9{v) = -v^- (D-57) 

1 - e-^— 

This is a lowered isothermal (Maxwell) distribution which vanishes at the escape velocity. 
For the fundamental mode, we then have 

/(-;;, t) = Ae-I^l* f e'^^ - e-f^^) , (D.58) 



where ^ is a normalization constant. Since the evaporation rate A is small, the system 
is in a quasi-stationary state given by the fundamental eigenfunction (King model). The 
total distribution slowly changes in amplitude (without change in form) as stars leave the 
system. Note that an extension of the King model to the case of fermions (or for the 
Lynden-Bell |23] theory of violent relaxation) has been obtained by Chavanis |40| . 

(i) The vortex model f38^ : in a cosmogonical context, it has been proposed [4T1|42] that 
large-scale vortices may have spontaneously emerged in the protoplanetary nebula, and 
have captured and accumulated dust particles, leading ultimately to planetesimals and 
planets. Chavanis [38] has used a Fokker-Planck approach to estimate the evaporation of 
dust from the vortices due to turbulence (with the final result that evaporation is negligible 
for relevant particles' sizes). This model is based on a Fokker-Planck equation of the form 
(j2.3p (where v plays the role of the position x) with a constant diffusion coefficient D in 
d = 2 dimensions. In that case, the function x{v) is given by 

pP4- r 2 

^[v) = ——l e-^'^wdw, (D.59) 
with x(0) = 0- This is easily integrated in 







..2 



Xiv) = ^ ( - 1 ) . (D.60) 
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leading to 



I3D Jq w 

This can be rewritten in the form of a series as 



X(v) = 4^ I - ( e^"^ - 1 ) dw. (D.61) 



' 2I3D ^ n\n ^ ' 

^ n=l 



This can also be written in the form 



X(^) = 1^ -j-\n(f3^]}, (D.63) 



where 



X gi 



Ei{x)= / -dt (D.64) 

J —oo 

is the exponential integral and 7 is the Euler constant. We see that 

for V large. The asymptotic evaporation rate ()D.48P is then 

X{V) ~ -I)/3V2e-^^'/2. (D.66) 

E Simpler expressions of the function G{v,u) 

Introducing a spherical system of coordinates, the function G{v, u) defined by (12. lip can 
be written 

(E.67) 

Using the identities 

"(sing)'^"^cig= ^^ffl , (E.68) 



r(i) 

j\~^<^°^\sir,ef~'i de = ' r (^^^^ V^/|_^(x), (e.69) 
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we obtain 



ci-2 



2 ) \Ru\v\ ) l-n2 2-1 Vl-^^^ 

In particular, we get 



G(Rw, n = r - - — -— —5- — -^e ^d-.^) ^ 

We also recall that 

Let us now consider particular dimensions of space d = 1,2,3 where the 
(|E.70P can be further simplified. 
• In d = 3, using the identity 

h/2{x) = a/— sinh(x), 
or directly integrating (IE.67p . we obtain 



This can also be written 



0(..„) 1 ' I ^ 



2(27r)3/2 Ru\v\ V 1 - u2 
In particular. 



g 2(l-u2) g 2(l-u2) 



2(27r)3/2 R^u \ 1- V? 
In d = 2, we obtain 



g 2(l+u) g 2(l-u) 



1 B f l3Ru\v\ 

1^ _ on _„2) r / ^ II 



G(u,n) = — 2(1-^2) 



27r 1 - u2 ^ V 1 - ^2 

In d = 1 , using the identity 

I-i/2{x) = W — cosh(a;) 
V vrx 



30 



or directly integrating ()E.67p . we obtain 

1 



G{v,u) = 
This can also be written 
G{v,u) = 

In particular, 

G{Ruj,u) -- 



(27r)i/2 V 1 - -u 



B _£(fl!i£!+H!) /BRulvl 
e 2(1-"^) cosh ' ' ' 



1-^2 



2(27r)V2 V l-'u2 



/5 



g 2(l-u^) _j_ g 2(l-u^) 



/3 



/3fl^(l-n) 



_£Hf(l+tO 



g 2(l+u) _|_ g 2(l-u) 



2(27r)i/2 V 1 - u2 

Finally, we note that for R — )■ +oo, the eigenfunction (12.12P is given by 



G{v,u) 



G{v,0) 

F Recurrence relations for the moments 

Let us introduce the moments of the distribution function 



du 
u 



Mn{t) 



f{v,t)\v\'^''dv, nGN 



(E.79) 



(E.80) 



(E.81) 



(E.82) 



and integrate Eq. (j2.3p with D = 1 against |?;|2",n > 1. After easy computations, one 
gets 

M^{t) = / Vf -nda- 2/3nM„(t) + 2n(2n + d- 2)M„_i(t). 



Doing the same for n = 0, we obtain 



= / Vf-nda, 



so that the foregoing expression can be rewritten 

M^(t) = R^'^M'Q{t) - 2(3nMn{t) + 2n(2n + d- 2)M„_i(t), n > 1. 
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